Mon. Not. R. Astron. Soc. 000, [TJg|(2010) Printed 28 October 2010 (MN M£X style file v2.2) 



Adaptive image ray-tracing for astrophysical simulations 

E. R. Parkin* 

Institut d'Astrophysique et de Geophysique, Universite de Liege, 17, Allee du 6 Aout, B5c, B-4-000 Sart Tilman, Belgium 
Accepted 2010 October 20. Received 2010 October 20; in original form 2010 July 20 



ABSTRACT 

A technique is presented for producing synthetic images from numerical simulations 
whereby the image resolution is adapted around prominent features. In so doing, 
adaptive image ray-tracing (AIR) improves the efficiency of a calculation by focusing 
computational effort where it is needed most. The results of test calculations show 
that a factor of > 4 speed-up, and a commensurate reduction in the number of pixels 
required in the final image, can be achieved compared to an equivalent calculation 
with a fixed resolution image. 
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1 INTRODUCTION 

Numerical models, such as hydro-dynamical simulations, 
have become a popular tool in astrophysical research. 
A common approach to extracting observable quantities 
from a numerical simulation is to trace the path of a 
ray through the simulation domain and sum-up the value 
of some quantity at discrete intervals. Typically this is 
done to integrate the column density through the sim- 
ulation domain or to solve the equation of radiative 
transfer to determine the emergent intensity of emission. 
Synthetic images produced via ray-tracing can be very 
useful for constraining model parameters, and examples 
■ of this can be found in a wide range of astrophysi- 
cal s tu dies, e.g. symbiot i c recu rrent nova (|Orlando et al.l 
120091: iDrake fc Orlandol |2010[). active galactic nuclei 
( Briiggen et al l 120091; iFalceta-Goncalves et al.l l2010l ) . star 



an image the resolution must be sufficiently high to ensure 
that the smallest scales of the simulation are well sampled. 
Yet there may be regions of a simulation which do not war- 
rant such a high level of sampling. Furthermore, a ray which 
passes through a highly refined region of the simulation do- 
main may not necessarily hold much useful information once 
it exits the grid. For instance, regions of high intrinsic emis- 
sion may be heavily absorbed such that the emergent inten- 
sity is negligible. With these details in mind, as well as the 
fact that as the computational requirements of simulations 
rise there will be an associated increase in the memory re- 
quired by datasets, more efficient approaches to extracting 
observable quantities are warranted. 



formation (iKurosawa et al. 20*041: Krumholz et al. 20071: 



Parkin et all 120091; lOffner fc KrumhohJ 120091: iPeters et al 



2010; Dougla s et al 



.201(1 ) 
20071; 



, jets (ISutherland fc Bicknelll 
Saxton et al.l |2010| ). starburst 



20071 ; iBonito et a". . 

galacti c" winds (| Cooper et alj 20081), and colliding winds b i 



naries (|Pittard fc Doughertv^OOd ; IPittard fe Parkinll2010l ). 

However, one drawback with using a fixed resolution image 
(i.e. same pixel size at all points) is that a large amount of 
computational effort can be expended on essentially blank 
regions as features of interest rarely fill an entire image. 

For example, in modern hydrodynamic simulations 
the grid/particle resolution is often adapted to features 
in the flow. In grid-based codes this is achieved using 
adap tive-mesh refinement (AMR - e.g. iBerger fc Oligerl 



19891). whereas sm oothed-particle hydrodynamics (SPH - see 
Monaghanl 1 19921 ) is inherently adaptive. When producing 
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This letter describes a method for adapting the reso- 
lution of a ray-traced image to the feature(s) of interest. 
In so doing the efficiency of a calculation is significantly 
improved. Adaptive image ray-tracing (AIR) takes advan- 
tage of the varying simulation resolution and magnitude of 
the extracted information encountered by a ray. In essence, 
the advantages that adaptive mesh refinement (AMR) has 
brought to grid-based hydrodynamical simulations are incor- 
porated into ray-tracing. The process is in a sense the reverse 
of su per-sampling (e.g. IWhittedl Il980l ; iGenetti fc Gordon! 
Il993l ). where instead of averaging over multiple rays, a single 
ray is subdivided. Test calculations show that, compared to 
a fixed resolution image, AIR provides an appreciable speed- 
up and a reduced number of pixels for the resulting image. 
The remainder of this letter is structured as follows: in § [5] 
the AIR technique is outlined, §[3] presents results from test 
calculations, and conclusions are presented in §|4] 
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2 ADAPTIVE IMAGE RAY-TRACING 

The basic principle behind producing a ray-traced image is 
to first discretize the plane of the sky into a uniform array 
of pixels and then follow the path of a ray for each respec- 
tive pixel. The AIR technique builds on this by taking an 
initially low resolution image and then adapting (increas- 
ing) the resolution around sufficiently prominent features of 
interest. This leads to a far more efficient calculation as com- 
putational effort is concentrated on producing an image of 
a desired variable(s). The structure of the AIR scheme is as 
follows: 

(i) Construct base image: After reading in the simu- 
lation file the base image resolution can be determined and 
the initial image constructed. As a guideline, the base image 
resolution can be set to sample the lowest resolution regions 
of the simulation domain. For example, for an AMR simula- 
tion setting the base image resolution equivalent to the base 
grid of the simulation is adequate. 

(ii) Ray-trace: Extract the desired information from the 
simulation domain by following the path of rays for each 
respective pixel. 

(iii) Scan the image to check if refinement is re- 
quired: This is a relatively straightforward process of cal- 
culating the truncation error for a given pixel. The step in 
resolution between adjacent image pixels is prevented from 
being more than one refinement level - this ensures that the 
edges of features are well sampled. 

(iv) Refine pixels: The pixels selected by the truncation 
error check should now be refined. The ij indexing of pixels 
cannot be preserved and therefore a heirarchical tree struc- 
ture is required. Hence, during this step the book-keeping 
must be performed and information about the neighbours, 
parent, and children of a pixel updated. 

(v) Loop: Repeat steps (ii)|(iv) until features in the im- 
age have been captured to the desired resolution and no 
more pixel refinement is required. For instance, the image 
resolution need not exceed that of the simulation. 

(vi) Integrate quantities and output: Once the final 
image has been constructed integrated quantities can be de- 
termined. An example of this would be broadband images 
in a given frequency /energy range. The sharpness of edges 
in the image can be improved by incorporating a super- 
resolution algorithm to int erpolate between adjacent pixels 
of differing resolution (e.g. IChu et al.|[2009h ■ Following this, 
the final step is to output the image. 

As a note of caution, when calculating an integrated 
spectrum from an image it is necessary to store the spec- 
trum for each individual pixel until the calculation is fin- 
ished. This can lead to a prohibitive overhead in memory 
requirements. However, this problem can be alleviated by 
writing the spectrum for each pixel to temporary files dur- 
ing the calculation. 

Considering that the majority of the computational ef- 
fort is expended on the ray-tracing step, effective speed- 
up can be achieved by parallelizing the AIR scheme on 
this step alone. This would involve distributing the list of 
pixels across processors during each ray-tracing sweep and 
then gathering the information for the refinement check. 
Alternatively, it is conceivable that such a code could 
be manufactured by taking an existing fixed-image ray- 



tracing code and inter facing it with a paralle l AMR library, 
e.g. PARAMESH jMacNeice et al.l |2000T), CHOMBO 



dColella et all I2009T), DAGH (|Parashar fc Brownel 1 19951 ). 



SAMRAI |Wissink et"aT1l200lh . For instance, AMR libraries 
typically have in-built functionality for parallel grid manage- 
ment, e.g.. refinement, cell list book-keeping, interprocessor 
communication, and load balancing. Therefore, ray-tracing 
could be handled on a pixel-by-pixel basis using an existing 
fixed-image ray-tracing code, and the parallel computation 
and grid management could be dealt with by functions avail- 
able in the AMR library. The process would be akin to the 
initial refinement sweep performed in a grid-based hydrody- 
namics code, with the difference that instead of populating 
a refined cell using the simulation initial conditions a re- 
fined pixel is populated using the results from a ray-tracing 
calculation. 



3 AN EXAMPLE APPLICATION 

To demonstrate the advantages of using AIR compared to 
ray-tracing with a fixed resolution image, calculations have 
been performed for a test case. For this purpose a ring of 
hot gas residing at the centre of a three-dimensional box has 
been simulated. The ring is aligned with the j/z-plane, has a 
constant thickness w = 1 x 10 17 cm, and has inner and outer 
radii of r\ = 3.5 x 10 17 cm and r = 4 x 10 17 cm, respectively. 
The box has dimensions x = y = z — ±5 x 10 17 cm. The 
gas density (g cm" 3 ), 



P = 



1 x 10" 25 



g cm 



: otherwise 



\x\ < w/2 



(1) 

where [3 = 2.5 x 10 9 g cm 1 and r yz — y 2 + z 2 . The gas 
temperature is 10 s K and 10 4 K inside and outside of the 
ring, respectively. 

The ring and the box are modelled o n an AMR grid con- 
structed using the FLASH code v3.1.1 (|Frvxell et al.ll2000l ; 
iDubev et aI1l2009h . which oper ates with the PARAM ESH 
block-structured AMR package |MacNeice et aLlbOOol ). The 
coarse grid consists of 4 3 blocks containing 8 3 cells. Re- 
finement is performed on density and temperature using 
4 nested grid levels such that the effective resolution is 
512 3 cells. 

For the AIR calculations the base image has a resolu- 
tion of 32 2 pixels and 4 nested levels of image refinement 
were used, giving an effective image resolution of 512 2 pix- 
els (equivalent to a fixed image calculation) . The refinement 
check is performed on the integrated 1-10 keV intrinsic X-ray 
flux using a mod i fied s econd-derivative interpolation error 
estimate (|Lohnenll987h . This is essentially a second-order 
central difference normalized by the sum of first-order for- 
ward and rearward differences, which in one dimension on a 
uniform mesh ifl, 
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\p i+ i - 2pi +Pi-i\ 
\p i+1 — Pi | + \pi -Pi-i\ 



(2) 



In iLohnerl 's formulation of the error estimator there is an ad- 
ditional term in the denominator added as a filter to prevent the 
refinement of ripples in hydrodynamic simulations. This term is 
not necessary for our purposes and, therefore, is not included. 
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where £ is the truncation error, p is the value of the pa- 
rameter on which image refinement is desired in pixel i (e.g. 
X-ray flux). The multidimensional generalization of Eq. [2] is 
found by taking all cross derivatives, which leads to, 



E[(l 



dp I , I 



dp 



-1/2 



Aar, 



r 

/ 



(3) 



where u and v are the image coordinates with indices i and 
j, respectively, and Ax a is the separation of nodes in the co- 
ordinate direction a. The partial derivatives are determined 
at pixel centres and the sums are carried out over coordinate 
directions. If £^ ^ £ C rit then the pixel is marked for refine- 
ment. Determining the optimal value for £ cr it for a given 
application requires some experimentation; for the test case 
values of £ cr it > 0.01 are effective at refining the ring whilst 
maintaining efficiency. However, further tests performed on 
images which contain more structure and varying degrees of 
contr ast revea l that £ cr it ~ 0.5 is a more appropriate start 
point. lL~ohner] 's error estimate is useful because it is bounded 
(i.e. ^ £ < 1) so that a preset refinement tolerance can 
be employed. It is also dimensionless, meaning that more 
than one image parameter can be used to check for refiment 
without encountering dimensioning problems. 

To calculate the intrinsic X-ray emission we assume so- 
lar abundances and use emissivities for optically thin gas 
in collisional ionization equilibrium obtained from look-up 
tables calculated from the MEKAL plasma code (jKaastral 
1 19921 ; iMewe et al.ll 19951 ) . 

The test problem places the observer viewing the ring 
face-on (parallel to the z-axis) at a distance of 1 kpc and 
traces rays through the hydrodynamic grid to calculate the 
intrinsic 1-10 keV flux. For comparison, calculations have 
been performed with a fixed image and with AIR using 
£crit = 0.2, 0.5, and 0.8 (Fig. [T] and Table E). The fixed 
image calculation is constrained by the fact that image pix- 
els must be small enough to sample the highest refined re- 
gions of the hydrodynamic grid, hence the resolution must 
be 512 x 512 pixels. In contrast, the AIR calculation begins 
with a base image of 32 x 32 pixels then locates the feature 
of interest (the ring in this case) and increases the image 
resolution respectively. This leads to far fewer pixels being 
used and thus a more efficient calculation. For example, the 
Ccrit = 0.8 AIR calculation takes ~ 1/7 the time and re- 
quires ~ 1/12 the image pixels, with a negligible error in 
the integrated flux of 0.01%. 

As a approximate rule of thumb, AIR will reduce 
the calculation time and pixel consumption by a factor of 
roughly the inverse of the image filling factor, e.g. in the test 
calculation the ring fills ~ 1/8 of the image. 



4 CONCLUSIONS 

A method has been presented for adaptively increasing the 
resolution of a ray-traced image around prominent features 
of interest. An initially low resolution image is scanned 
and relevant pixels are refined to increase the image resolu- 
tion. The results of test calculations show that considerable 
speed-up (a factor of ~ 4 — 7), and a commensurate reduc- 
tion in the number of pixels required for the final image (a 



Table 1. Adaptive image ray-tracing test calculation results. £ cr ; t 
is the critical error estimate used to flag pixels for refinement (see 
Eq|3j, Fx is the intrinsic 1-10 keV flux, t is the time taken for the 
calculation. The error is calculated as (F Xe xact - ^XrayV-Fxexact, 
where -Fxexact = 1.399586 X 10 — 19 erg s _1 cm~ 2 is the intrinsic 
flux summed from the hydrodynamic grid and Fxray is the ray- 
traced value for a given calculation. 



Calculation 


£crit 


Fractional Error 


t 


Pixels 








(s) 




Fixed 




~ 0. 


9604 


262144 


AIR 


0.2 


2 X 10~ 6 


2627 


40512 


AIR 


0.5 


1 X 10~ 4 


1374 


23384 


AIR 


0.8 


1 x 10~ 4 


1339 


22748 





Figure 1. Results from the gcrit = 0.8 test calculation showing 
intrinsic 1-10 keV X-ray flux (top panel) and the image pixel 
mesh (lower panel). The images show a spatial extent of u = v = 
±5 X 10 17 cm - large tick marks correspond to 1 X 10 17 cm. Note 
that for the adopted viewing angles (u,v) = (y, z). For further 
details of the test calculations see §[3] and Table [T] 



factor of ~ 6 — 11), can be achieved compared to an equiva- 
lent calculation with a fixed resolution image. In conclusion, 
adaptive image ray-tracing (AIR) improves the efficiency of 
a calculation by focusing computational effort on extracting 
desired information. 
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